home *** CD-ROM | disk | FTP | other *** search
/ Complete Linux / Complete Linux.iso / xwindows / demos / xfract_1.z / xfract_1 / xfractint-1.06 / fpu087.c < prev    next >
C/C++ Source or Header  |  1992-09-28  |  4KB  |  186 lines
  1. /* Fpu087.c
  2.  * This file contains routines to replace fpu087.asm.
  3.  *
  4.  * This file Copyright 1991 Ken Shirriff.  It may be used according to the
  5.  * fractint license conditions, blah blah blah.
  6.  */
  7.  
  8. #include <fractint.h>
  9. #include <mpmath.h>
  10. #include <math.h>
  11.  
  12. extern int DivideOverflow;
  13.  
  14. /*
  15.  *----------------------------------------------------------------------
  16.  *
  17.  * xxx086 --
  18.  *
  19.  *    Simulate integer math routines using floating point.
  20.  *    This will of course slow things down, so this should all be
  21.  *    changed at some point.
  22.  *
  23.  *----------------------------------------------------------------------
  24.  */
  25.  
  26. void FPUaptan387(double *y, double *x, double *atan)
  27. {
  28.     *atan = atan2(*y,*x);
  29. }
  30.  
  31. void FPUcplxmul(struct complex *x, struct complex *y, struct complex *z)
  32. {
  33.     double tx;
  34.     tx = x->x * y->x - x->y * y->y;
  35.     z->y = x->x * y->y + x->y * y->x;
  36.     z->x = tx;
  37. }
  38.  
  39. void FPUcplxdiv(struct complex *x, struct complex *y, struct complex *z)
  40. {
  41.     double mod,tx,ty,yxmod,yymod;
  42.     mod = y->x * y->x + y->y * y->y;
  43.     if (mod==0) {
  44.     DivideOverflow++;
  45.     }
  46.     yxmod = y->x/mod;
  47.     yymod = - y->y/mod;
  48.     tx = x->x * yxmod - x->y * yymod;
  49.     z->y = x->x * yymod + x->y * yxmod;
  50.     z->x = tx;
  51. }
  52.  
  53. void FPUsincos(double *Angle, double *Sin, double *Cos)
  54. {
  55.     *Sin = sin(*Angle);
  56.     *Cos = cos(*Angle);
  57. }
  58.  
  59. void FPUsinhcosh(double *Angle, double *Sinh, double *Cosh)
  60. {
  61.     *Sinh = sinh(*Angle);
  62.     *Cosh = cosh(*Angle);
  63. }
  64.  
  65. void FPUcplxlog(struct complex *x, struct complex *z)
  66. {
  67.     double mod,zx,zy;
  68.     mod = sqrt(x->x*x->x + x->y*x->y);
  69.     zx = log(mod);
  70.     zy = atan2(x->y,x->x);
  71.  
  72.     z->x = zx;
  73.     z->y = zy;
  74. }
  75.  
  76. void FPUcplxexp387(struct complex *x, struct complex *z)
  77. {
  78.     double pow,y;
  79.     y = x->y;
  80.     pow = exp(x->x);
  81.     z->x = pow*cos(y);
  82.     z->y = pow*sin(y);
  83. }
  84.  
  85. /* Integer Routines */
  86. void SinCos086(long x, long *sinx, long *cosx)
  87. {
  88.     double a;
  89.     a = x/(double)(1<<16);
  90.     *sinx = (long) (sin(a)*(double)(1<<16));
  91.     *cosx = (long) (cos(a)*(double)(1<<16));
  92. }
  93.  
  94. void SinhCosh086(long x, long *sinx, long *cosx)
  95. {
  96.     double a;
  97.     a = x/(double)(1<<16);
  98.     *sinx = (long) (sinh(a)*(double)(1<<16));
  99.     *cosx = (long) (cosh(a)*(double)(1<<16));
  100. }
  101.  
  102. long Exp086(x)
  103. long x;
  104. {
  105.     return (long) (exp((double)x/(double)(1<<16))*(double)(1<<16));
  106. }
  107.  
  108. #define em2float(l) (*(float *)&(l))
  109. #define float2em(f) (*(long *)&(f))
  110.  
  111. /*
  112.  * Input is a 16 bit offset number.  Output is shifted by Fudge.
  113.  */
  114. unsigned long ExpFudged(x, Fudge)
  115. long x;
  116. int Fudge;
  117. {
  118.     return (long) (exp((double)x/(double)(1<<16))*(double)(1<<Fudge));
  119. }
  120.  
  121. /* This multiplies two e/m numbers and returns an e/m number. */
  122. long r16Mul(x,y)
  123. long x,y;
  124. {
  125.     float f;
  126.     f = em2float(x)*em2float(y);
  127.     return float2em(f);
  128. }
  129.  
  130. /* This takes an exp/mant number and returns a shift-16 number */
  131. long LogFloat14(x)
  132. unsigned long x;
  133. {
  134.     return log((double)em2float(x))*(1<<16);
  135. }
  136.  
  137. /* This divides two e/m numbers and returns an e/m number. */
  138. long RegDivFloat(x,y)
  139. long x,y;
  140. {
  141.     float f;
  142.     f = em2float(x)/em2float(y);
  143.     return float2em(f);
  144. }
  145.  
  146. /*
  147.  * This routine on the IBM converts shifted integer x,FudgeFact to 
  148.  * the 4 byte number: exp,mant,mant,mant
  149.  * Instead of using exp/mant format, we'll just use floats.
  150.  * Note: If sizeof(float) != sizeof(long), we're hosed.
  151.  */
  152. long RegFg2Float(x,FudgeFact)
  153. long x;
  154. int FudgeFact;
  155. {
  156.     float f;
  157.     long l;
  158.     f = x/(float)(1<<FudgeFact);
  159.     l = float2em(f);
  160.     return l;
  161. }
  162.  
  163. /*
  164.  * This converts em to shifted integer format.
  165.  */
  166. long RegFloat2Fg(x,Fudge)
  167. long x;
  168. int Fudge;
  169. {
  170.     return em2float(x)*(float)(1<<Fudge);
  171. }
  172.  
  173. long RegSftFloat(x, Shift)
  174. long x;
  175. int Shift;
  176. {
  177.     float f;
  178.     f = em2float(x);
  179.     if (Shift>0) {
  180.     f *= (1<<Shift);
  181.     } else {
  182.     f /= (1<<Shift);
  183.     }
  184.     return float2em(f);
  185. }
  186.